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We develop and test Quantum Monte Carlo algorithms which use a "twist" or a phase in the 
wave function for fermions in periodic boundary conditions. For metallic systems, averaging over 
the twist results in faster convergence to the thermodynamic limit than periodic boundary conditions 
for properties involving the kinetic energy with the same computational complexity. We determine 
exponents for the rate of convergence to the thermodynamic limit for the components of the energy 
of coulomb systems. We show results with twist averaged variational Monte Carlo on free particles, 
the Stoner model and the electron gas using Hartrec-Fock, Slater- Jastrow, three-body and backflow 
wavefunction. We also discuss the use of twist averaging in the grand canonical ensemble, and 
numerical methods to accomplish the twist averaging. 
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Almost all quantum Monte Carlo (QMC)calculations 
in periodic boundary conditions have assumed that phase 
of the wavefunction returns to the same value if a parti- 
cle goes around the periodic boundaries and returns to its 
original position. However, with these boundary condi- 
tions, delocalized fermion systems converge slowly to the 
thermodynamic limit because of shell effects in the fill- 
ing of single particle states. In this paper we explore an 
alternative boundary condition: one can allow particles 
to pick up a phase when they wrap around the periodic 
boundaries. 



can be restricted to be in the range: 



*(ri + Lx, r 2 , • • •) = e <9 -*(ri,ra, • ■ •)• 
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The boundary condition 9 = is called periodic bound- 
ary conditions (PBC), 9 = 7r anti-periodic boundary con- 
ditions (ABC) and the general condition with 9^0, 
twisted boundary conditions (TBC)a. 

In periodic boundary conditions, the Hamiltonian is 
invariant with respect to translating any particle around 
the periodic boundaries. According to Bloch's theorem, 
this implies that any solution can be characterized by 
a given twist angle. The twist angle also has a physi- 
cal origin: consider a toroidal geometry. One can either 
rotate the toruso and go into rotating coordinates, or 
add a magnetic fluxcl to the center of the torus. The 
physical properties will be unchanged. In both cases one 
can transform away the perturbation by applying TBC 
with the twist angle given by 9 = mR 2 ui/h for rotation 
and 9 = e<f>/(ch) for magnetic flux. A torus is topologi- 
cally equivalent to periodic boundary conditions, so that 
a non-zero twist will be allowed in periodic boundaries. 
The twist is a degree of freedom, or boundary condition, 
that can be varied to enable a finite system to approach 
the thermodynamic limit more quickly or to make de- 
tailed studies of the properties of the quantum state. 

If the periodic boundaries are used in all three direc- 
tions, each dimension can have an independent twista 
Hence, in 3D, the twist is a three component vector, 
9i with i = {1,2,3}. The free energy and hence all 
equilibrium properties are (triply) periodica in the twist: 
F(9i + 27r) = F(9i) so that each component of the twist 



7T < 



< 7T. 



(2) 



For systems with a real potential (e.g. no magnetic field), 
one can further restrict the twist to be in the range [0, ir]. 

For a degenerate Fermi liquid, finite-size shell effects 
are much reduced if the twist angle is averaged over. We 
call this twist averaged boundary conditions (TABC)lI. 
This is particularly important in computing properties 
that are sensitive to the single particle energies such 
as the kinetic energy and the magnetic susceptibility. 
By reducing shell effects, much more accurate estima- 
tions of the thermodynamic limit of these properties can 
be made. What makes this even more important is 
that the most accurate quantum methods have compu- 
tational demands which increase rapidly with the num- 
ber of fermions. Examples of such methods are exact 
diagonalizationB (exponential increase in CPU time with 
N), variational Monte Carld3(VMC) with_wavefunctions 
having backflow and three-body termaiffl (increases as 
./V 4 ), and transient-estimate and released-node Diffusion 
Monte Carlo methodso (exponential increase with N). 
Methods which can extrapolate more rapidly to the ther- 
modynamic limit are crucial in obtaining high accuracy. 
Twist averaging is especially advantageous for stochas- 
tic methods (i.e. QMC) because the averaging does not 
necessarily slow down the evaluation of averages, except 
for the necessity of doing complex rather than real arith- 
metic. 

The use of twisted boundary conditions is common- 
place for the solution of the band structure problem for 
a periodic solid. Band structure methods begin by as- 
suming the wavefunction factors into single particle or- 
bitals characterized by a lattice momentum. Then in 
order to calculate properties of an infinite periodic solid, 
properties must be averaged-by integrating over the first 
Brillouin zone. Baldereschiliil pointed out that in an in- 
sulator, in integrating over the Brilliouin zone, one can 
with high accuracy replace the integral with a "special 
k-point." This was generalized to a grid of k-pointsEj. 
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Twisted boundary conditions has been-discussed in con- 
nection with polarization of insulatorsa; we do not con- 
sider that here. The use of twisted boundara-jconditionsiii 
common in the analysis of lattice modelsBo. GammelMI 
showed using perturbation arguments for certain lattice 
models why it will converge faster to the thermodynamic 
limit—and applied it to calculating optical properties. 
Grosta studied size effects in the Hubbard model with 
exact diagonalization and showed TABC gives exact re- 
sults in the grand canonical ensemble for non-interacting 
systems. 

Though twisted boundary conditions have a long his- 
tory within quantum physics, their use in quantum 
Monte Carlo (QMC) has been limited. In continuum 
QMC, RajagopaO used special k-points to reduce fi- 
nite size errors for calculations of insulators within VMC. 
Their use in Diffusion Monte Carlo (DMC) results were 
restricted to PBC and ABG-in order to work with real 
wavefunctions. Kralik et aZ£3 used generalized bound- 
ary conditions to compute the momentum distribution 
with VMC. for silicon (a semiconductor) and Filippi and 
CeperleyOj for lithium (a metal). This was done in order 
to enlarge the number of momentum vectors to probe the 
behavior near the fermi surface. The use of TABC within 
QMC has not been further developed. 

The applications considered here are for systems with 
a fermi surface, i.e. metals. We begin by discussing 
the method for non-interacting (NI) fermions. Fermi 
liquid theory asserts that the spectrum of states of in- 
teracting liquids are intimately related to those of the 
non-interacting fermion systems, hence a detailed anal- 
ysis for the NI system carries over to strongly interact- 
ing fermi liquids. We then discuss interacting systems 
in the Hartree-Fock approximation: the electron gas and 
the Stoner model. In the Stoner model, we show how 
TABC can be used to determine a polarization phase 
transition at non-zero temperature. Results for TABC 
are given for the interacting electron gas using a pair 
product and backflow wavefunction in 3d. The electron 
gas system has been previously treated with an extrap- 
olation method based on Fermi liquid theory. We show 
that TABC gives the same results in the thermodynamic 
limit and verify the applicability of the NI analysis, in 
particular to examine how the energy depends on the 
twist of a given system size. We then present VMC re- 
sults of the polarization energy of the electron gas using 
the new method. In future publications we will study 
the low density properties of the electron gas using this 
technique. The Appendix discusses some details arising 
in the implementation of TABC in QMC. 



I. NON-INTERACTING FERMIONS 

In a non-interacting homogenous system with PBC, 
the single particle states are plane waves: exp(ikr)?7(<7) 
where r\ is a spin function. For simplicity, we always as- 



sume the simulation cell is a cube (or square in 2D) of 
side L. To satisfy the twisted boundary conditions, the 
wave vectors obey: 

k n = (2™ + 9)/L (3) 

where n is an integer vector. These states have energy 
E a = (h 2 /2m)\i n . The ground state in the canonical en- 
semble consists of the N lowest energy states; the many- 
body wavefunction is a determinant of those states. In 
this section we will ignore spin, since for a non-interacting 
system, spin modifies the results only by doubling the de- 
generacy of each level. Figure [l] shows the occupation of 
states for 13 spin-less fermions in 2D for 9 = i. e. with 
periodic boundary conditions, and also with a non-zero 
twist. The occupied states lie within a circle centered at 
the origin with radius w kp = 2{itp) 1 / 2 . 
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FIG. 1. Momentum distibution for 13 spinless fermions in 
a 2D square with side L — 2n. The top panel shows the oc- 
cupied states (closed symbols) and empty states (open sym- 
bols) with zero twist (circles, PBC) and a twist equal to 
27r(0.3, 0.15) (triangles). The circle shows the infinite system 
fermi surface. The bottom panel shows the occupied states 
with TABC. The colored regions show the occupied region 
for the lowest level (middle square), the third level, up to the 
outermost 13*' 1 level. 

Figure ^| shows the relative error in energy versus the 
number of fermions with PBC. The energy converges 
slowly to the exact result. One sees "cusps" in the curve 
at certain values of N. These occur at closed-shell values 
of N, e.g. the state depicted in Fig. [j] for PBC is a closed 
shell since states related by symmetry are either all filled 
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or all empty. For large N the curve is "quasi-random" , 
with an envelope decaying algebraically as N~ v . 

We find numerically that the exponent of the decay of 
the relative error of the energy is approximately v = 1.33 
in 2D and v = 1 in 3D (see Table I) . To characterize the 
approach to the thermodynamic limit, we introduce two 
different measures. Defining the relative scaled error 5n 
as: 

E N =E 0O {l + N- v 5 N ), (4) 

we define a = max \6n\, b — (Sn) and c = {(6n — b) 2 ) 1 / 2 . 
Table I shows estimates of these coefficients and ex- 
ponents obtained numerically by examining values of 
10 < N < 10 4 . 




FIG. 2. Relative error of the energy versus number of par- 
ticles with PBC (A) and TABC (□) in 2D and 3D. The points 
shown are only those where the relative error has a local max- 
imum. Curves are shown only for N < 100. 

Now consider twisting the boundary conditions, i. e. 
using a nonzero phase. This displaces the set of k- vectors 
as shown in the top panel of Fig. |l|. Aside from a set of 
special twists having zero measure, the energy levels will 
no longer be degenerate. (When we sort the states to 
decide the filling, all states will have a different energy.) 
This is because inversion symmetry and rotational sym- 
metry through 90° is broken. This breaking of symme- 
tries and absence of degeneracies is a crucial difference 
between TBC and PBC. 

At critical values of the twist, when a filled and empty 
state have the same energy, the occupation of the states 
changes. The condition for the degeneracy is that k lie 



on a plane bisecting and perpendicular to the line joining 
the origin with an integer vestor; precisely the Laue con- 
dition for the Bragg planesEi of the reciprocal lattice of 
the supercell. In fig. || is shown the dependence of total 
energy on the twist angle for a fixed number of particles. 
One sees cusps as the filled states cross the Bragg planes. 
The dependence is similar to the band energy of a simple 
metal. Later, we will discuss this band structure for an 
interacting system. The bandwidth Ebw (the spread of 
energy values fig. ||) is defined as: 

E 2 BW = (2ir)- d J d6{E{6) - £ M ) 2 (5) 

depends on the number of particles and scales as Ebw oc 
N~ v where the exponent is the same as describes the 
convergence of the KE in PBC. 

There are several alternative procedures by which the 
twist angle can be varied: i) one can average the twist 
over all possible values, ii) the twist can become a dy- 
namical variable and iii) special values of the twist could 
be used. Of these approaches, none are right or wrong 
in general; which method approaches the thermody- 
namic limit faster depends on the order of the phase 
in question, whether fermi liquid, ferromagnetic or anti- 
ferromagnetic. However, to compute a variety of prop- 
erties for a metallic systems, we find the TABC best re- 
duces size effects. 
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TABLE I. Coefficients of the asymptotic decay of the error 
in the relative NI energy, v is the exponent of the decay. The 
exponents have been determined from the numerical data and 
are accurate to about 0.02. The amplitude was determined 
numerically by examining the values for 10 < N < 10000. 
The coefficients are defined as a = max(|5jv|), b — (5n), 
c = ([8m — b] 2 ) 1 ^ 2 . n is the number of phase angles used for the 
summation in each dimension: 0i = 2ni/n for i — 1, . . . ,n. 
n — 1 corresponds to PBC; d is the dimensionality. O is the 
property: T, the kinetic energy; V, the Hartree-Fock poten- 
tial energy of the electron gas; N, the number of particles in 
the TA-GCE method. 
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II. TWIST AVERAGING 



The twist average of a property A is defined by: 

(A) = (27r)- d r d9(^Rj)\A\^(Rj)) (6) 



where it is assumed that the wavefunction i/j(R, 9) is nor- 
malized for each 9. The momentum distribution, n(k) is 
a key property to calculate for delocalized quantum sys- 
tems. A discontinuity in n(k) at the Fermi surface is 
responsible for the validity of Fermi Liquid Theory for 
metals. The kinetic energy is the second moment of the 
momentum distribution: 



T N = {h 2 /2m) J dkk 2 n{k). 



(7) 



Let us analyze the momentum distribution for NI 
fcrmions in the canonical ensemble. For any given twist 
9, the N lowest energy states from those given by Eq. 
|§| are occupied. But any value of k can only be reached 
by a unique combination of (n, 9) if 9 is restricted by 
Eq. |^. This proves that the averaged momentum dis- 
tribution is a constant for states which can be reached 
by some combination of (n, 9) and zero otherwise. Hence 
within TABC, the set of filled states comprises a "solid 
volume" bounded by a Fermi surface. In contrast, for a 
single twist value, the momentum distribution is a point 
set. The total volume in k-space inside the Fermi surface 
is precisely (27r) d p, just as it is in the thermodynamic 
limit, so the constant is determined by the normalization 
condition: 



/ 



dkn(k) = 1. 



(8) 
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As mentioned above, the Fermi surface is a subset of 
the Bragg planes. For N particles the occupied-states 
comprise the union of the first N Brillouin zoncsEj. The 
(N + l) th electron will go in the (N + l) th zone, an area 
formed by planes surrounding the N th zone. Figure |l| 
shows the momentum distribution of 13 spinless fcrmions 
in 2D using TABC. 

In ID, TABC gives the exact momentum distribution 
because the normalization condition determines every- 
thing. In higher dimensions the fermi surface is not per- 
fectly circular (spherical) as shown in fig. [I]. However, 
one can see that n(k) is much closer to a disk than the 
momentum distribution obtained with PBC. A perfect 
fermi surface (no finite size corrections) in any dimen- 
sion, can be obtained by allowing variations in the parti- 
cle number and working in the grand canonical ensemble 
as we discuss below. 



FIG. 3. Dependence of the energy of NI unpolarized 
fermions on the twist angle for JV=54 in 3D. The solid line 
shows the energy along the (100) direction, dotted line along 
the (110) direction and dashed line along the (111) direction. 
The curves are piecewise quadratic, with a cusp when the oc- 
cupation of the states changes. We refer to the r.m.s. spread 
of energy values as the bandwidth. 



Figure || shows the convergence of the error of the ki- 
netic energy within TABC versus the number of parti- 
cles. Numerical estimates of the relative error are given 
in Table I. One sees a dramatic improvement in the con- 
vergence with respect to PBC. The exponents governing 
the decay rate are larger and the errors are up to two 
orders of magnitude smaller for system sizes in the range 
N ks 100. Note that the TABC kinetic energy must ap- 
proach the exact energy from above. This is because the 
shape of a given volume with the smallest moment of in- 
ertia is a sphere, so that the distorted shape shown in fig. 
]l] has a higher energy. Also shown in Table I is depen- 
dence of the error on the number of twist values in the 
average. One needs from 16 to 32 values of 9 along each 
axis to achieve the full reduction in size effects (better 
than a percent accuracy in the relative error of the size 
effects). In the Appendix are discussed the relative mer- 
its of performing the average on a grid versus sampling 
the twist values from a uniform distribution. 



Let us now examine how the potential energy converges 
with PBC and TABC. This will give us some idea of how 
two particle correlations are affected by the boundary 
conditions since the potential energy is a particular inte- 
gral over the pair correlation function. The calculation 
performed below is particularly simple for a power law 
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potential, v(r) = r~". In particular, we examine the po- 
tential energy of an electron gas {n = 1) computed using 
the NI wavefunction (Hartree-Fock approximation) . The 
NI trial function is valid for high density when the kinetic 
energy dominates the potential energy. The potential en- 
ergy (using the Ewald image potential) is conveniently 
evaluated in Fourier space as a sum over the structure 
factor: 



V 



Np 



$^«kOSk-i) + i\r«jif 



(9) 



TABC. For all values of N and twists the potential en- 
ergy of the finite systems approaches that of the infinite 
systems from below. Twist averaging serves to make the 
decay more regular but does not reduce its overall mag- 
nitude which is determined by a charge interacting with 
the correlation hole in the nearby supercell. Similar ef- 
fects are expected for other two particle quantities. The 
smoother convergence obtained with TABC should allow 
for more accurate extrapolation to N — > oo. 



where vm is the Madelung energy of a charge interacting 
with itself and «k is the Fourier transform of the interpar- 
ticle potential. For a 1/r potential, v% = 27r(<i — l)/k x . 
The values of k in the sum are k = 27rn/L where n is 
an integer vector. For the NI wavefunction, the structure 
factor at wave- vector q is proportional to the probability 
that after we have displaced a filled state by q we are 
still in a filled state. 



q)> 



(10) 



k<k' 



where the sum is over occupied states and the average is 
over twisted boundary conditions. 
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FIG. 4. Relative error in the evaluation of the potential 
energy for an electron gas using the Hartree-Fock wavefunc- 
tion for N spinless electrons. The solid line shows TABC, the 
dashed line PBC. 

Shown in fig. |i| is the convergence of the potential 
energy versus the number of particles using PBC and 



III. GRAND CANONICAL ENSEMBLE 

The use of TABC in the grand canonical ensemble 
(GCE) gives the exact single particle occupations for NI 
particles, as shown by GrostLa within the Hubbard model. 
Suppose (a, N) is a label of the quantum states, both for 
the number of particles and for other quantum numbers 
such as the momentum and let E a ^{6) be the energy of 
this state. Then the probability of a given state in the 
GCE is proportional to exp(—(3(E a! N{&) — Nfj,)) where fi 
is the chemical potential. In the ground state, /3 — > oo, 
the occupied many-body state will be the one minimizing 
E a ,N{6) ~ Nfi. Thus N can depend on 9. 

We now show that for a NI system, TABC within the 
GCE gives exact single particle properties; i.e. there are 
no finite size effects. Suppose the single particle energy 
levels are et- Then the probability of occupying the N 
states ei, . . . ,e„ is exp(— 5Zfc=i fti e k{9) — /•*])■ In the oc- 
cupation number njt basis, this probability distribution 
factorizes as: 



nf 



n k e 



+ (l-n k ) 



(11) 



so the probability of state k being occupied is precisely 
the fermi distribution law n k = (exp(/3(efc — fx)) + 1) . 
As the twist angle is varied over its range, each momen- 
tum state of the infinite system occurs precisely once. 
Hence the averaged occupation number is precisely what 
it would be in the thermodynamic limit. The distortion 
of the fermi surfaces observed in the lower panel of fig. |l| 
is a consequence of using the canonical ensemble. 

The momentum distribution and hence the kinetic en- 
ergy will be exactly equal to their infinite system values. 
Other properties may have finite size corrections, only 
the single particle properties are guaranteed to be exact. 
We call this procedure the twist average grand canonical 
ensemble (TA-GCE). 

With this procedure one does not have a fixed number 
of particles since for a given twist and fermi wave vector, 
the number of occupied states will vary. The fluctuations 
in the number of particles is closely related to a famous 
problem in analytic number theory, "Gauss' circle prob- 
lem" , to determine the number of lattice points inside a 
circle of area A as its radius tends to infinitjO. As Gauss 
posed the problem, the center of the circle was fixed on 
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a lattice site while in TA-GCE the circle is placed rap* 
domly (the shifted circle problem). It has been shownE3 
that in 2D one obtains the following fluctuation in the 
particle number for TA-GCE: 

lim ((N - {Lkpf/Air) 2 ) 1 ' 2 <x N 1 / 4 . (12) 

N—*oo 

Our numerical estimates for the convergence are shown 
in Table I. Note that the Gauss circle problem differs 
from the convergence of the energy in the canonical en- 
semble since the energy is the average second moment of 
the lowest N points, not the number of points in a circle. 

The preceding discussion refers to NI systems, how- 
ever, fermi liquid theory asserts that the low lying ex- 
cited states of an interacting system are in one-to-one 
correspondence withJie NI states. Hence, as argued by 
GrossEl and GammelEl, TA-GCE is likely to reduce finite 
size effects substantially for interacting systems. 

For an interacting many-body system, one difficulty in 
using the GCE is the need to optimize the wavefunction 
at each twist value; in particular to pick out which or- 
bitals should be occupied. For an isotropic system having 
a spherical fermi surface, the order of filling the states is 
simple. For a metal with a non-spherical Fermi surface, 
the usual procedure is to determine the filling of single 
particle states according to a mean-field theory such as 
the Kohn-Sham method in density functional theory. If 
one uses the same procedure within QMC, the Fermi sur- 
face will be substantially unchanged. Hence, it would be 
better to try other fillings, choosing the one which mini- 
mizes En — Nfx. 

There are other problems in using the GCE for charged 
systems in periodic boundary conditions. Usually the 
positive compensating charge, either a uniform back- 
ground or a fixed array of charged nuclei, is at a fixed 
density. But if the number of electrons fluctuates, the 
periodic cell can have a net charge, which causes prob- 
lems in calculating the long-range potential. Although 
very promising, we do not consider the TA-GCE method 
further in this paper. The following examples are for the 
canonical ensemble. 



IV. THE STONER MODEL 

To test the utility of TABC for determining a phase 
transition, we simulatjcd the Stoner model, an analyti- 
cally soluable modeOE-3 for strongly correlated-fermions, 
particularly related to itinerant magnetismEj. The 
Stoner model differs from NI fermions by the addition 
of a contact repulsive potential: Yli<j 9^( r ij)- The en ~ 
ergy is evaluated within the mean field (Hartree-Fock) 
approximation using the NI wavefunction. For a spin 
(1/2) system, the potential energy is E = Enj + gn^ri^. 
In 3D the energy at zero temperature in the thermody- 
namic limit is: 



E/N = A^^- [(1 + C) 5/3 + (1 - C) 5/3 + Kl - C 2 )] 

(13) 

where b = ^(3vr 2 p)- 2 / 3 and A = h 2 /2m. For b < 1.111 
the system has an unpolarizcd ground state and for 
b > 1.3228 the ground state is ferromagnetic. For in- 
termediate couplings, the ground state has a partial spin 
polarization, at zero temperature, similar to the behavior 
suggestecOeZI for the electron gas at low density. 

We performed a MC simulation of the 3D NI sys- 
tem, within the occupation representation. The state 
variables of the system consist of occupation numbers 
(both spin and wavevector) and the twist angles. We 
make Monte Carlo moves consisting of flipping spins and 
moving the spatial occupation while keeping the parti- 
cle number fixed. If the new state is not already occu- 
pied, the move was accepted with probability equal to 
exp(-/3(e ne w(0) - e o id{0))). Fig. | shows the conver- 
gence of magnetization distribution versus the number of 
fermions. We found that one could determine the phase 
transition within a few percent accuracy in the coupling 
constant b using only 54 fermions. Such accuracy within 
PBC would require thousands of fermions because of the 
strong shell effects. 
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FIG. 5. Magnetization as a function of b for the 3D Stoner 
model. Circles are the exact results, diamonds for 54 electrons 
with PBC, squares for 54 electrons with TABC and triangles 
for N=200 with TABC. The temperature was T = 0.224Sf 
and a 8 3 grid of twist values was used. 

V. ELECTRON GAS 

We now present results for the 3-D electron gas as 
a test of TABC on a correlated many-body continuum 
system. The electron gas is a very important model 
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in condensed matter physics being the basis for the 
density functional theory method of electronic structure 
computations .tl The phase diagram of the electron gas 
at low density is still not resolvecO. The wavefunctions 
we use are the most accurate known for a continuum 
many-body system with optimized 2-body, 3-body and 
back flow terms. Variational MC (VMC) and Diffusion 
(DMC) are QMC methods appropriate to zero tempera- 
ture, and Path Integral (PIMC) to T > 0. In this paper, 
we will only discuss VMC and DMC. PIMC will be dis- 
cussed in future publications. 

In VMC, one assumes an analytic form for a trial 
function ^t(R) where R symbolizes the 3iV coordi- 
nates. Then one samples |\1/t(-R)| 2 using a random 
walku. An upper bound estimate to the exact ground 
state energy is the average of the local energy El(R) = 
^t(R)~ 1 'H^t(R) over the random walk. An accurate 
trial wavefunction is obtained from the NI wavefunction 
by multiplying by pair correlation terms. The pair prod- 
uct or Slater-Jastrow wave function is: 

* 2 (i?) = Det(e l ^)e~^«J U{r ' j) . (14) 




FIG. 6. Energy versus the number of electrons for the 3D 
electron gas using SJ wavefunction. The circles and connect- 
ing line are TABC with 10 3 twist values. The filled A's are 
using PBC. The arrow shows the extrapolation to an infinite 
system made using the PBC calculations using Eq. 

One can determine the correlation factor u(r) either 
with an analytic argument or by a minimizing of the en- 
ergy and/or the variance of an assumed form. We use a 
parameter-free analytic form so that the systematic size 
and twist effects are not masked by noise in the trial 
function itself. For the electron gas very accurate ana- 
lytic form for uLr) based on the energy within the RPA 
approximatiortJ has as low an energy as those with opt 
timized parameters. We used optimized Ewald sumsEHI 
both for the potential and for the correlation factor so as 
to have the correct long wavelength behavior. 



Fig. |6| shows the convergence of the VMC energy ver- 
sus the number of particles within TABC and PBC. One 
can see the convergence to the thermodynamic limit that 
was found within the NI system is also evident within 
VMC using the SJ wavefunction. In addition to the 
Slater-Jastrow trial functions, we also have, used opti- 
mized backflow-three body functions(BF3B)El which give 
a more accurate description of the low-density electron 
gas. (They pick up about 3/4 of the remaining corre- 
lation energyu and break certain symmetries of the SJ 
wavef unct ions . ) 

In the Diffusion Monte Carlo (DMC) method, one 
starts with a trial function and uses exp(— tH) to project 
out the ground state using a branching random walk. 
Fermi statistics pose a significant problem to the pro- 
jection method, since exact methods such as transient 
estimate or release-node QMC suffer a exponential loss 
of efficiency for large numbers of particles. For this rea- 
son the approximate fixed-node method is normally used. 
Using the fixed-node method, one obtains the best upper 
bound to the energy consistent with an assumed sign of 
the wave function. Both the FN method and the exact 
transient estimate caiL-be generalized to treat complex- 
valued trial functions.!^! These methods are called the 
fixed-phase and released phase QMC. We have tried both 
of these approaches using TABOlj 

VI. FERMI LIQUID THEORY AND TABC 

Fermi liquid theory (FLT) for metallic systems allows 
both a method to extrapolate to the thermodynamic 
limit and a way of understanding the twist dependence 
of the QMC results. According to Landau, the low-lying 
excitations of an interacting system are in close relation 
to those of the NI system. 

E = E + J dkfe k e(k) + J dkdk'SrikSnvfCk, k') + . . . 

(15) 

where Sn^ is the deviation of the quasi-particle occupa- 
tion from the ground state, and e(k) and /(k, k') are 
1 and 2 quasiparticle energy functionals. For simplicity 
we have not indicated dependence on spin. The energy 
functionals are usually further expanded about the fermi 
surface in spherical harmonics and applied to calculate 
properties in the thermodynamic limit. Here we discuss 
how to apply FLT when the "excitation" is caused by 
the boundary conditions. For example, we consider the 
momentum distribution shown in fig. |l|. The change in 
the energy caused by the non-circular shape should also 
be given by Eq. [l5| 

We can analyze this dependence by comparing the en- 
ergies of the non-interacting system within a given twist 
with the interacting system. This is done in Fig. |t]. One 
sees a linear relation between the two energies, confirm- 
ing that for these wavefunctions, Fermi liquid theory is an 
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appropriate description. The slope is the effective mass 
of the quasi-particles. Then since the NI infinite energy 
is known this gives a way of determining the interacting 
system energy. 

Previous calculations on the electron gas have usee 
other application of FLT, the extrapolation methodc 
In this method, one calculates accurate ground state en- 
ergies for a sequence of particle numbers, and determines 
the effective mass, potential correction and infinite sys- 
tem energy by fitting these energies to the relation: 



E N = £00 + (m/m*)AT N + eN~ 



(16) 
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where AT/v is the deviation of the NI kinetic energy from 
the infinite system and v is the exponent for the potential 
energy given in Table I. Figure || shows that the estimate 
of the infinite system energy obtained using TABC and 
the extrapolation method agree within errors. This is 
reassuring, but expected, since both are based on FLT. 

However, applied in practice, there are significant ad- 
vantages to the TABC method. The fitting method re- 
quires well converged runs for at least three different val- 
ues of N. This can be difficult, for example when the 
unit cell is large. For example, suppose one is doing a 
simulation of bcc hydrogen where one is limited to val- 
ues of N equal to 2K 3 = {2,16,54,128,250...}. The 
last 3 values are reasonable to simulate, however, if the 
unit cell contained 10 electrons, for example, the extrap- 
olation, would be very difficult without further approx- 
imations. There is also the problem that convergence 
and trial functions for large N may differ from that for 
small N. For example, one typically determines m/m* 
and e within VMC and then applies the extrapolation 
using the DMC energies. There are other problems with 
extrapolation. For example, suppose we have a partially 
spin polarized system with say ri\ spin up particles and 
n,2 down particles. Using PBC there are no necessarily 
larger, closed shell systems with precisely the same ratio 
£ = (m - n 2 )/{n 1 + n 2 ). 

Within TABC, there is the possibility of obtaining re- 
sults in the thermodynamic limit, using only quantities 
computed for a single value of N. Figures [7||| demon- 
strate that the 
Other methodsE! 



dnetic contribution is well controlled, 
are needed to make the potential cor- 
rection also for a fixed N. Use of TABC should make 
this more accurate as well, as demonstrated in fig. f|. 



FIG. 7. Plot of the VMC energy versus the deviation of 
the NI energy from Eoo. Each phase angle is plotted sepa- 
rately. Simulations are for N = 54 unpolarized 3D electron 
gas at r s = 50. Since the excitations are linearly related, 
Fermi liquid theory describes the phase angle dependence. 
Upper points are the Slater-Jastrow wavefunction, lower ones 
are done with the backflow-3-body trial function. Effective 
masses (the slope) are respectively 1 and 0.61 for the two 
trial functions. 
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FIG. 8. Energy versus spin polarization for the 3D elec- 
tron gas at r-s=40 and N = 54. The dashed line is with the 
Slater- Jastow function, the solid line with backflow-3-body 
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Calculations are done using VMC with 10 
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VII. OTHER PROCEDURES USING TWISTED 
BOUNDARY CONDITIONS 

Dynamical twist method. As an alternative to 
TABC, one can let the twist angle be a dynamical vari- 
able (DTBC) and be determined self consistently. By 
dynamical is meant that the probability of a given twist 
is proportional to its free energy: 

P{9) cx exp(-/3F(0)) = J2 e M-f3E a (6)) (17) 

a 

where E a (9) is the energy of the state a with twist 9. 
This distribution of twist angles could be attained within 
QMC by enforcing detailed balance on moves of the twist 
angle during a random walk. The expectation of an op- 
erator A will be: 

(A) D tbc = \\ dee-e F W(y a (6)\A(6)\y a (6)}. (18) 

At zero temperature, a special set of twist angles, those 
which have the lowest energy, will be singled out, so that 
the ground state energy will be Edtbc = vc&a${E(ff)). 
Clearly Edtbc < E^, in contrast to the TABC which 
gives an upper bound. In general, the dynamical energy, 
Edt converges as slowly to the thermodynamic limit as 
does the PBC energy: the exponents and fluctuations 
are the same. The dynamical twist method is not an 
improvement over PBC for approaching the thermody- 
namic limit for a metallic system at temperatures much 
lower than the fermi temperatures. 

Although the dynamic twist method is not satisfactory 
for a Fermi liquid at low temperatures, for certain lattice 
models such as an antiferromagnetic Heisenberg model, 
it can be a definite improvement over PBC and TABC. 
This is because one can only establish a defect-free Neel 
ordering if the unit cell is commensurate with the bound- 
ary conditions. For the Heisenberg model on a triangular 
lattice, the ground state has a given twist per lattice spac- 
ing. Using DTBC allows one to establish this twist value 
automatically, without imposing it in advance. This is 
equivalent to the classical variable cell method where the 
dimensions and aspect ratio of the supercell of a crystal 
become dynamical variables, so that one-.can determine 
the most stable crystal lattice structured instead of ex- 
amining each crystal structure explicitly. 

Special points For each value of N there exists a set of 
twist values for which E(9) — E^. One can determine 
these special twist values for the NI system and then per- 
form simulations at only one of those twist values for the 
interacting system, thereby getting rid of single particle 
size effects. This is similar to the special k-point method 
of Baldereschitll for insulators where a single k-point is 
determined by symmetry, thus allowing one to replace 
an integral over the Brillouin zone with evaluation at a 
single k-point. 31is method was used within QMC by 
Rajagopal et al.tZL 



The special k-point method is not appropriate for a 
metal because of the discontinuity in properties at cer- 
tain twist values. One cannot replace the average by 
a single point because the Fermi surface is not given in 
advance by symmetry, and can change between the inter- 
acting and non-interacting wave functions. In addition, 
it is not expected that the same twist values appropriate 
for the NI energy, will be appropriate for other quantities 
such as the potential energy, or spin susceptibility. It is 
better to have a method which can give a spectrum of 
properties correctly, rather than only the kinetic energy. 
As we discuss in the Appendix, TABC does not impose 
an excessive computational burden, and is to be preferred 
over using only special twist values. 



VIII. CONCLUSIONS 

Note that there is a significant difference in the ef- 
ficiency for stochastic (QMC) methods versus explicit 
methods such as density functional theory or exact di- 
agonalization in regards to the TABC. In explicit meth- 
ods, computations for each twist require an equal amount 
of computer time so that averaging over Nq twist values 
will take roughly Nq times as long as a single twist value. 
One can use inversion and rotational symmetry to reduce 
this, so that for a grid of 16 3 points, TABC will require 
only 165 twists. On the other hand, in QMC all the twists 
reduce the statistical error of the average. The twists 
are simply three more degrees of freedom on top of the 
3iV coordinate variables to be averaged over. However, 
one must also take into account startup costs associated 
with each twist, such as re-optimizing the trial wavefunc- 
tion or equilibrating the random walk. Neglecting these 
startup costs, there is no loss of statistical efficiency in 
performing TABC so that the gain in reducing the sys- 
tematic error is free. This is examined in more detail in 
the Appendix. 

There are many examples where twist averaging can ef- 
fect considerable improvement over the use of PBC. We 
have been able to perform quite accurate calculations pL 
the polarization energy of the 2 and 3D electron gasEBl 
and of liquid 3 He using backflow wavefunctions with on 
the_prder of 100 fermions. As pointed out by Ortiz ct 
al.EZl calculations on such small systems in PBC have con- 
siderable systematic errors. One property that could be 
computed more accurately with TABC is the estimation 
of fermi liquid parameters by calculating particle-hole 
excitatior£3. TABC can reduce the shell effects which 
caused much difficulty in that calculation. A related ex- 
ample is in rGomputation of the charge response of the 
electron gastll where a considerable effort was made to 
cancel out effects of the PBC. We apa presently study- 
ing the electron gas confined to a slabEa to determine the 
work function and surface energy of a metallic surface. 
The filled states consist of a set of disks, each of which 
will have a certain occupation number. By doing twist 
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averaging we have shown reduced size effects with respect 
to PBC. 

Experimental systems are at a non-zero temperature. 
For NI systems one occupies the states with probability 
given by the Fermi-Dirac distribution. Because the Fermi 
function at non-zero temperature is a continuous func- 
tion, the convergence to the thermodynamic limit will be 
much faster, even with PBC. However in practice, one is 
interested in electronic systems close to the ground state; 
the relevant quantity is the thermal deBroglie wavelength 
of the electron: %/ '{m e k b T) 1/2 « 32A at T = 300K. Be- 
cause this length is usually larger than the simulation cell 
in QMC, the localization of the density matrix does not 
help at reducing fermion finite size at these temperatures. 
Using a non-zero temperature just to achieve faster con- 
vergence to the thermodynamic limit is not practically 
useful. However, TABC is extendable to finite temper-, 
ature PIMC simulations using the fixed-phase methodtil 
and will reduce size effects at low temperature. 

The TABC method is likely to be valuable for all QMC 
calculations in systems with a fermi surface. The calcula- 
tions on the electron gas demonstrate that even though 
it may be a little slower per step of the random walk, 
it is better to do TABC than a larger system with peri- 
odic boundary conditions because TABC converges much 
faster to the thermodynamic limit. The overall efficiency 
of any numerical method is ultimately judged by the com- 
puter time needed to reduce systematic and statistical 
error below a given value. TABC is effective in reduc- 
ing the systematic errors and thus improve the overall 
efficiency. 
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tional resources were provided by the NCSA. We thank 
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APPENDIX 

Here we discuss numerical details of implementing 
twist averaged boundary conditions. Many of the 
changes caused by twisted boundary conditions arise 
from the need to have complex wavefunctions. Although 
the wavefunction and energies in special cases are real 
(e.g. PBC and ABC), complex functions are needed for 
general twist angles. There is a factor of roughly 2.5 in 
additional CPU time to do the arithmetic to evaluate 
the determinant and its derivatives. The actual impact 
on the total speed is smaller than this because the cal- 
culations of two-particle quantities such as the potential 
energy and correlation factors are still done with real 
arithmetic; the actual penalty of working with non-zero 
twist depends on the number of particles and the type 
of the trial wavefunction. However, as we have discussed 



earlier, even a factor of 2.5 in computer time is worth- 
while if one is able to approach the thermodynamic limit 
quicker, since QMC methods scale as N v with 1 < v < 4 
or, in the case of exact fermion methods, as exp(7iV). If 
TABC saves going to larger N, the additional time doing 
complex arithmetic is well justified. 

There are several alternatives for performing the twist 
averaging: 

1. Evaluate as WiEi using a grid defined by points 
6i with weights yu (with W{ = 1) in the region 
specified by Eq. ||. 

2. Sample the twist during the QMC random walk 
and take the average. 

3. A combination of the two approaches: working on a 
grid that is augmented with random displacements. 

As we discuss below, all three methods are satisfactory; 
there is no fundamental difference in efficiency. The 
choice of whether to sample or use a grid is primarily 
based on convenience and programming considerations 
and only secondly on efficiency. We note that all meth- 
ods are easy to parallelize. 

Grid averaging. First, we must address the question of 
which grid and integration rule to use. Since all prop- 
erties are periodic with respect to the twist, the grid 
should be a Bravais lattice with equally weighted points. 
One must keep in mind that the properties, though pe- 
riodic and continuous, have discontinuous derivatives at 
the Bragg planes. Unless grid points can be located on 
these planes (which is difficult to achieve in practice), the 
integration error will go as eg oc A9 2 oc N g 2 ^ D where Nq 
is the number of grid points. Numerically, we find that 
16 3 grids are needed for an accuracy of 1CU 3 (see table 
I). This slowly convergent, systematic error is the main 
drawback of the grid integration method. For an insu- 
lator with a large enough gap to excitations, properties 
would be analytic for all 9 since the occupation of single 
particle states will not change as a function of twist angle 
and the grid error would converge exponentially fast. 

Once a grid is chosen, one can use symmetry (e.g. in- 
version and rotation through 90°) to reduce the num- 
ber of grid points and give them a weight (wi with 
w i = 1) proportional to their multiplicity. It 
is easy to show that the optimal amount of computer 
time at each grid point should be chosen proportional to 
Wi/^Oi) 1 / 2 where the MC efficiency at 0i is defined as 
C(0j) = l/(var(E(9i))cputime). We have found on the 
calculations of the electron gas described earlier, that 
this efficiency is independent of the twist angle except at 
the special PBC and ABC points where real functions 
can be used. Even though we have symmetry and can 
integrate over a reduced set of twist values, we must in- 
tegrate longer at high multiplicity points since they con- 
tribute more to the average. Hence, the symmetry does 
not significantly reduce the needed amount of CPU time. 
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Since the calculations at different twist angles are un- 
correlated, one can easily show that the efficiency of cal- 
culating the twist averaged energy is given by the rela- 
tion: 

r 1/2 =X>c 1/2 (^). (19) 

Hence, the overall efficiency of the TABC energy is an 
average of the efficiencies of the individual twist calcula- 
tions and is higher than that of the slowest converging 
twist angle. The additional averaging over twist angle 
costs nothing in efficiency. 

However, this discussion did not take into considera- 
tion start-up costs at each twist angle, such as the need 
to reoptimizc the trial wave function at a new twist value, 
and equilibration costs. By equilibration, we mean that 
whenever the twist angle is changed, enough random walk 
steps must be taken so that the configurations are sam- 
pled from the new twist value. During this equilibra- 
tion, averages cannot be taken. These computational 
costs cause a decrease of efficiency by the factor (useful 
time) / (total time) and are the main extra computational 
penalty of the TABC method within QMC. Since the 
startup time will scale with the number of needed grid 
points, using the above estimate of the systematic er- 
ror, ec, we find that the needed startup time scales as 
e G^ 2 w hile the time to achieve equivalent statistical er- 
ror scales as (T 2 . Hence, for very precise calculations 
(e — » 0) and d < 4, startup costs can be neglected. 

Twist sampling. Now consider the second alternative, 
where the twist angle is sampled during the random walk. 
With this method, we do not have to decide on a grid in 
advance and there is no systematic error of a finite grid. 
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